Proximity-induced sub-gaps in Andreev billiards 
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We examine the density of states of an Andreev billiard and show that any billiard with a finite 
upper cut-off in the path length distribution P{s) will possess an energy gap on the scale of the 
Thouless energy. An exact quantum mechanical calculation for different Andreev billiards gives 
good agreement with the semiclassical predictions when the energy dependent phase shift for An- 
dreev reflections is properly taken into account. Based on this new semiclassical Bohr-Sommerfeld 
approximation of the density of states, we derive a simple formula for the energy gap. We show that 
the energy gap, in units of Thouless energy, may exceed the value predicted earlier from random 
matrix theory for chaotic billiards. 
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Studies of sub-gap transport in hybrid superconduetors 
are an important starting point for the design and simu- 
lation of superconducting nanoscale devices. For a given 
inhomogeneous structure, a fundamental question is the 
existence or otherwise of a finite density of states at the 
Fermi energy. This is important experimentally, since the 
sub-gap spectrum determines the tunneling conductance 
of an N-S contact. The abihty to address this question is 
also important theoretically, since a well-posed problem 
of this kind provides a testing ground for complementary 
(and occasionally competing) theoretical techniques. 

One important class of structures, for which a general 
analysis might be forthcoming, are known as Andreev 
billiards. These are formed when a classically-chaotic 
normal dot is placed in contact with a superconductor 
[0, H, |[ ^, |[ ^, ^ . Consider a ballistic two dimensional 
normal dot of area A, with the mean level spacing of the 
isolated normal system 5 — at the Fermi energy 

Ep. If a superconductor of width W and bulk order 
parameter A is placed in contact with such a billiard 
(see Fig. p, then the question of interest is whether or 
not an energy gap on the scale of the Thouless energy 
Et = ^ 1^, exists in the sub-gap spectrum E < A, 
when S << Et < A. The number of open channels in 
the S region is the integer part of M = and the 

energy levels of the Andreev billiards are the positive 
eigenvalues E (measured from the Fermi energy) of the 
Bogoliubov-de Gennes equation [p^. 
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FIG. 1: Geometry. 



An initial study of this problem based on random ma- 
trix theory^, concluded that a classically chaotic billiard 
possesses an energy gap, whereas an integrable system is 
gapless. As a consequence it was suggested that the exis- 
tence of such a gap could be used to distinguish between 
integrable and chaotic systems. Later it was shown that 
billiards with classically mixed phase space do possess 
a smaller gap comparing to the fully chaotic systems 
Studying the existence of the gap for chaotic billiards 
Ihra et al. have drawn the attention to the role of 
the non-diagonal terms of the scattering matrix which 
could explain the disagreement between the random ma- 
trix theory and the Bohr-Sommerfeld semiclassical ap- 
proximation. The non-universal feature of the excitation 
spectra of the Andreev billiards has been recently studied 
by Ihra and Richter in Ref. 



The aim of this Letter is to identify for the first time, a 
pseudo-integrable billiard which does possess an energy 
gap. Moreover, this gap can be larger than that pre- 
dicted for chaotic billiards on the scale of Et- To this 
end, we study the quasi-particle spectrum of the ballis- 
tic structure shown in Fig. |l|, for different values of the 
angle a and length d. Our key result is illustrated in 
Fig. 1^, which shows that for d — 0, the counting func- 
tion or integrated density of states N{E) — n{E') dE' 
(where n{E) is the density of states for the Andreev 
billiard) exhibits an energy gap -Egap varying between 
Egap ~ 0.5Et and Eg^p ~ 1.5Et as a varies from 80° 
to 45° (for the details of the calculation see below). In 
contrast, as shown in Fig. ^, for d ^ 0, no such energy 
gap exists. To understand this result, we start from the 
formula relating the integrated density of states N{E) to 
the path length distribution P{s) (which is the classical 
probability that an electron entering the billiard at the 
N-S contact exits after a path length s) derived by using 
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FIG. 2: The integrated density of states, N{E) as a function 
of E (in units of A) for d = and angle a = 80°, 60°, 45°. 
For each angle a the stair type lines correspond to the exact 
diagonalization of the Bogoliubov-de Gennes equation, while 
the solid lines are obtained by using Eq. (|l]). The parameters 
are M = 55.5, A/Ef = 0.015. 




FIG. 3: The integrated density of states, N{E) as a function 
of E (in units of A) for d / and angle a = 80°, 60°, 45°. 
In each case, d varies such that the area of the normal dot is 
fixed at ^ = 5W^ for the angles a = 80° , 60° , 45° . This choice 
results in approximately the same number of states below A 
for different values of a. The meaning of the lines and the 
parameters are the same as those in Fig. ^ 



the Bohr-Sommerfeld approximation, 
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W -|j is the coherence length and 



the distribution P{s) is normahzed to one, /g P{s) ds 

density o\ 
y) yielding 



1. The density of states n{E) = can be found from 

Eq, " 
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This expression reduces to that of derived in a different 
way by Melsen et al. ^, Lodder and Nazarovj^, fhra 
et al. 0, in the hmit E << A and « A, but more 
generally incorporates an energy dependent path length 
correction pdj ]. 

To test the semiclassical expression, we perform an 
exact (numerical) diagonalization of the Bogoliubov-de 
Gennes equation by matching the wave functions at the 
N-S interface. This results in a secular equation includ- 
ing the scattering matrix, SpjE) of the normal billiard 
opened at the N-S interface 
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S,s{E)S:s{~E) =0, (4) 



where 



SME) = [Q{E) + K{E)D{E)]-' 
[Q*{E)^K{E)D{E)], 

D{E) = [l-So{E)][l + Sp{E)]- 



(5) 



Here Q and K are diagonal matrices with elements 

Qnm{E) = iSnrnQniE) and Knm{E) = 5nmkn{E), 

where qn{E) = kpy^ 



are the trans- 



verse wavenumbers of the electron in the S region and 
k„(E) = kp^/l + S — -T^T are the transverse wavenum- 

It is 



XE) = k,^l + ^-^ 
bers of the electron in the S region when A = 0. 
assumed that the Fermi wavenumber, kp — \j2mEp/h?' 
is the same in the S and N regions. All the matrices 
are M by M dimensional. Finally, Sq was calculated 
using Bessel functions in the wedge part of the normal 
region, and including evanescent modes. In a different 
context, the same type of normal billiard was studied by 
Kaplan and Heller A secular equation similar to 

Eq. (^ was derived by Beenakker for SNS systems 
but there instead of Scs{E) the scattering matrix So{E) 
of the normal billiard appears. In this sense, the matrix 
Scs{E) can be regarded as an effective scattering ma- 
trix. Note that in Andreev approximation p5[| qn ~ fcn, 
and one then finds Scs{E) = So{E). In contrast the sec- 
ular equation (H) for Andreev billiards is valid outside 
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Andreev approximation. Using the unitarity of Sq{E), 
one can show that the foUowing equation 
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[Q{E)+K{E)D*{-E)]-'}) = 0,(6) 



has the same zeros as those of Eq. (Q), and is a more 
suitable form for numerical calculations (the left hand 
side of Eq. (^) is a real function instead of a complex 
one). 

Figures § and |^ show both exact results and an eval- 
uation of the semiclassical expression of N(E) given by 
Eq. (|l|) for Andreev billiards shown in Fig. |] with differ- 
ent parameters One can see that the results of the 
two methods are in good agreement over the entire en- 
ergy range < A. It is important to note that inclusion 
of the energy dependence of the phase shift was essential 
to obtain such a good agreement. A small deviation can 
be seen only at energies close to the value of A and for 

Both the semiclassical and the exact calculations reveal 
the presence of an energy gap when d = 0, while no gap 
appears in case of d 7^ 0. The origin of this result can 
be traced to the existence of a finite upper cut-off Smax 
in the path length distribution P{s). From Eq. (^) it 
is obvious that if the path length of the electron has a 
maximum value, then the energy spectrum will possess a 
gap. The condition for determining -Egap is Sn{E) > Smax 
for n > 0. Expanding the arccos term in Eq. (|^) to first 
order in E/ A, one finds that the energy gap is given by 
the following simple equation: 
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where A is the area of the billiard and W is the width of 
the superconductor lead. 

As examples. Fig. | shows P(s) for a = 80°, 60°, 45° 
along with the corresponding semiclassical density of 
states n{E) given by Eq. (^). From the figure it is seen 
that the upper cut-off Smax of P(s) equals to 2W for 
a = 80°, 60°, and s^ax = 2V2W for a = 45°. It is in- 
teresting to note that for a = 45° (more generally for 
a = 90° /k with an integer k > 0) the density of states 
n{E) has a singularity at some energy E < A while for 
other values of a no pronounced singularities exist. In 
Bohr-Sommerfeld approximation, it can be shown that 
this singularity is related to the singularity of the path 
length distribution P(s). The classical trajectories of the 
electron resulting in a singularity for a — 45° are similar 
to those moving in a lead with length d+W and a = 90°. 
The nature of the singularities of the density of states in 
this rectangular shape of Andreev billiard was studied in 
Ref.|l|. The reason for the disappearance of these singu- 
larities for a 90° /k (k = 1,2,3,.. .) is consistent with 
the points made in this reference. 
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FIG. 4: The path length distribution P{s) for Andreev billiard 
of Fig. |l| when d = for a = 80°, 60°, 45°. Insets show the 
corresponding semiclassical density of states n[E) (in units of 
1/5) using Eq. (^). The same parameters was used as those 
in Fig. |. 



To check the expression (^ we have determined the 
energy gap from exact calculations for several angles a. 
For the bilhard of Fig. |l| with d = 0, one can show ana- 
lytically that the upper cut-off s„ 
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where k = 1, 2, 3, ... . Using Eqs. @} and (||) the gap is 
plotted as a function of the angle a in Fig. || together 
with the exact results. The numerical results (full circles 
in Fig. |5|) agree very well with the predictions of Eq. (Q). 
One can see that decreasing a the energy gap in units of 
Thouless energy increases and tends to a finite value as 
a ^ 0. From Eq. (||) and Eq. (||) it can be seen that Smax 
increases with decreasing a, whereas Egap/A decreases. 
The Thouless energy Et also decreases, since the mean 
level spacing 5 becomes smaller for larger area A of the 
normal region. The two effect together result in a finite 
value of E^^^p/Et at a = 0. 

In the limit ^0 << Smax (for example, in the semi- 
classical limit of large M), the expression (|^) of the gap 
reduces to 
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In Fig. H the dashed curve is the result from this formula. 
One can see that for all angles a this limiting result for 
the energy gap is larger than the finite M value. Note 
that the area is A = ^ cot a, and therefore, Egap/Er 
is only a function of a and independent of W when ^0 < < 
5max- It can also be shown that PgapZ-Ex 7r^/4 sa 2.47 
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density of states (even for energy levels close to the value 
of A) and for the energy gap. Finally we have shown 
that the energy gap on the scale of the Thouless energy 
can be much larger than the value O.GEt predicted from 
random matrix theory]^ for chaotic billiards. 
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FIG. 5: The gap -Egap (in units of Thouless energy) as a func- 
tion of angle a for Andreev billiards with d = 0. The solid and 
the dashed curves are the semiclassical results from Eq. (^) 
when AI = 55.5, and from Eq. (^ i.e. when M oo, respec- 
tively. The circles are the results of the exact diagonalization 
when M = 55.5. 



as a — + 0. Finally, it is worth mentioning that for Andreev 
billiard of Fig. yj when d = the energy gap can be much 
larger for small angle a than the value of O.QEt predicted 
by Melsen et al. [|| for the chaotic billiard. 

It is interesting to mention that besides the billiard 
studied in this work (see Fig. |^ for d = 0) there are 
a number of different shapes of normal dots in which 
P{s) exhibits finite upper cut-off resulting in a sub-gap. 
Examples for such dots are shown in Fig. ^ (provided 
a + P < 180° in the case of quadrangle). 
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FIG. 6: A triangle (a), and quadrangle (b) shape of normal 
billiard having a finite upper cut-off of the length distribution 
Pis). 

In conclusion, we have shown that a billiard with a 
finite upper cut-off in the path length distribution P{s) 
will possess an energy gap in the density of states on the 
scale of the Thouless energy. By including the energy 
dependent phase shift of the Andreev reflection, a new 
expression for the density of states has been given within 
the frame work of the semiclassical Bohr-Sommerfeld ap- 
proximation. We have also derived a formula for the 
energy gap. To check these results we have performed 
an exact diagonalization of the Bogoliubov-de Gennes 
equation for different Andreev billiards. The results of 
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